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Abstract. Boundary integral equations and Nystrom discretization provide a powerful tool 
for the solution of Laplace and Helmholtz boundary value problems. However, often a weakly- 
singular kernel arises, in which case specialized quadratures that modify the matrix entries 
near the diagonal are needed to reach a high accuracy. We describe the construction of four 
different quadratures which handle logarithmically-singular kernels. Only smooth boundaries 
are considered, but some of the techniques extend straightforwardly to the case of corners. Three 
are modifications of the global periodic trapezoid rule, due to Kapur-Rokhlin, to Alpert, and 
to Kress. The fourth is a modification to a quadrature based on Gauss-Legendre panels due to 
Kolm-Rokhlin; this formulation allows adaptivity. We compare in numerical experiments the 
convergence of the four schemes in various settings, including low- and high-frequency planar 
Helmholtz problems, and 3D axisymmetric Laplace problems. We also find striking differences 
in performance in an iterative setting. We summarize the relative advantages of the schemes. 



1. Introduction 

Linear elliptic boundary value problems (BVPs) where the partial differential equation has 
constant or piecewise-constant coefficients arise frequently in engineering, mathematics, and 
physics. For the Laplace equation, applications include electrostatics, heat and fluid flow, and 
probability; for the Helmholtz equation they include the scattering of waves in acoustics, elec- 
tromagnetics, optics, and quantum mechanics. Because the fundamental solution (free-space 
Green's function) is known, one may solve such problems using boundary integral equations 
(BIEs). In this approach, a BVP in two dimensions (2D) is converted via so-called jump rela- 
tions to an integral equation for an unknown function living on a ID curve [7]. The resulting 
reduced dimensionality and geometric simplicity allows for high-order accurate numerical solu- 
tions with much more efficiency than standard finite-difference or finite element discretizations 

The BIEs that arise in this setting often take the second-kind form 

(1.1) cr{x)+ [ k{x,x)a{x)dx' = f{x), xG[0,r], 







where [0,T] is an interval, where / is a given smooth T-periodic function, and where is a 
(doubly) T-periodic kernel function that is smooth away from the origin and has a logarithmic 



singularity as — )• x. In order to solve a BIE such as (1.1) numerically, it must be turned into 
a linear system with a finite number unknowns. This is most easily done via the Nystrom 
method [20l I18j. (There do exist other discretization methods such as Galerkin and collocation 
[18j : while their relative merits in higher dimensional settings are still debated, for curves in 
the plane there seems to be little to compete with Nystrom [71 Sec. 3.5].) However, since the 
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split into Lp, ijj explicit 


split into if, unknown 


global 


• Kresst yiTJ 


• Kapur-Rokhlin [T3] 


(periodic trapezoid rule) 




• Alpert Xl 






o QBX* [15] 


panel-based 

(Gauss-Legendre nodes) 


o Helsing [HI US] 


• Modified Gaussian (Kolm-Rokhlin) [16] 
o QBX* [H] 



Table 1. Classification of Nystrom quadrature schemes for logarithmically- 
singular kernels on smooth ID curves. Schemes tested in this work are marked 
by a solid bullet ("•"). Schemes are amenable to the FMM unless indicated with 
a f. Finally, * indicates that other analytic knowledge is required, namely a local 
expansion for the PDE. 



fundamental solution in 2D has a logarithmic singularity, generic integral operators of interest 
inherit this singularity at the diagonal, giving them (at most weakly-singular) kernels which we 
write in the standard "periodized-log" form 

(1.2) k{x, x') = (fix, x') log (^4 sin^ + ^/,(x, x') 

for some smooth, doubly T-periodic functions f and ip. In this note we focus on a variety of 
high-order quadrature schemes for the Nystrom solution of such ID integral equations. 

We first review the Nystrom method, and classify some quadrature schemes for weakly- 
singular kernels. 

1.1. Overview of Nystrom discretization. One starts with an underlying quadrature scheme 
on [0,T], defined by nodes {xi}^^ ordered by < rci < 2:2 < 3^3 < • • • < xn < T, and 
corresponding weig hts {u^ilili. This means that for g a smooth T-periodic function, 

rT N 



/ g{x)dx y]wi5f(xi) 
-'0 ,_i 



holds to high accuracy. More specifically, the error converges to zero to high order in N. Such 
quadratures fall into two popular types: either a global rule on [0, T], such as the periodic 
trapezoid rule [181 Sec. 12.1] (which has equally-spaced nodes and equal weights), or a panel- 
based (composite) rule which is the union of simple quadrature rules on disjoint intervals (or 
panels) which cover [0, T]. An example of the latter is composite Gauss-Legendre quadrature. 
The two types are shown in Figure [l] (a) and (b). Global rules may be excellent — for instance, 
if g is analytic in a neighborhood of the real axis, the periodic trapezoid rule has exponential 
convergence [IHl Thm. 12.6] — yet panel-based rules can be more useful in practice because they 
are very simple to make adaptive: one may split a panel into two smaller panels until a local 
convergence criterion is met. 



The Nystrom method for discretizing (1.1) constructs a linear system that relates a given 



data vector / = {/i}^^ where = f{xi) to an unknown solution vector cr = {ai}iLi where 
ai ^ <7{xi). Informally speaking, the idea is to use the nodes {xj}^^ as collocation points where 
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( 1.1 ) is enforced: 



(1.3) (y{xi)+ I k{xi,x')a{x') dx' = f{xi), z = 1, 

Jo 

Then matrix elements {aj,j}f^=i are constructed such that, for smooth T-periodic a, 

(1.4) / k{xi,x')a{x') dx' ^ 2,^1 j ^i^j) ■ 

Jo ,=1 



Combining (1.3) and (1.4) we obtain a square hnear system that relates a to /: 

N 

(1.5) CTi + '^aijaj = fi, i = l,...,N . 



node. We write (1.5) in matrix form as 



In a sum such as (1.5), it is convenient to think of Xj as the source node, and Xi as the target 



(1.6) (T + Aa = f. 

A high-order approximation to cr(x) for general x € [0, T] may then be constructed by interpo- 
lation through the values cr. 

If the kernel k is smooth, as occurs for the Laplace double-layer operator, then the matrix 
elements 

(1.7) aij = k{xi,Xj)wj 

lead to an error convergence rate that is provably the same order as the underlying quadrature 



scheme \18\ Sec. 12.2]. It is less obvious how to construct the matrix A = {ctij} such that (1.4) 



holds to high order accuracy in the case where k has a logarithmic singularity, as in (1.2). The 
purpose of this note is to describe and compare several techniques for this latter task. Note 
that it is the assumption that the solution a{x) is smooth (i.e. well approximated by high-order 
interpolation schemes) 

1.2. Types of singular quadrature schemes. We now overview the schemes presented in 
this work. It is desirable for a scheme for the singular kernel case to have almost all elements be 



given by (1.7), for the following reason. When N is large (greater than 10^, say), solving (1.6) 
via standard dense linear algebra starts to become impractical, since 0{N^) effort is needed. 
Iterative methods are preferred which converge to a solution using a small number of matrix- 
vector products; in the last couple of decades so-called fast algorithms have arisen to perform 
such a matrix- vector product involving a dense N x N matrix in only 0{N) or 0(A^log A^) time. 
The most well-known is probably the fast multipole method (FMM) of Rokhlin-Greengard [9j, 
but others exist [3l [51 [25] . They use potential theory to evaluate all N sums of the form 

N 

(1.8) '^k{xi,Xj)qj, i = l,...,N 



where the qj are interpreted as charge strengths. Choosing qj = Wjaj turns this into a fast 
algorithm to evaluate Act given cr, in the case where A is a standard Nystrom matrix (1.7). 
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(a) 




(c) 



** ** 
* 



Figure 1 . Example smooth planar curve discretized with = 90 points via (a) 
periodic trapezoid rule nodes and (b) panel-based rule (10-point Gauss-Legendre; 
the panel ends are shown by line segments). In both cases the parametrization 
is polar angle t £ [0,27r] and the curve is the radial function f{t) = 9/20 — 
(1/9) cos(5t). (c) Geometry for 2D Helmholtz numerical examples in section 7.2 



and 7.3, The curve is as in (a) and (b). Stars show source locations that generate 



the potential, while diamonds show testing locations. 



Definition 1.1. We say that a quadrature scheme is FMM- compatible provided that only 0{N) 



elements {flij} differ from the standard formula (1.7). 



An FMM-compatible scheme can easily be combined with any fast summation scheme for the 



sum (1.8) without compromising its asymptotic speed. Usually, for FMM-compatible schemes. 



the elements which differ from (1.7) will lie in a band about the diagonal; the width of the band 
depends only on the order of the scheme (not on A'"). All the schemes we discuss are FMM- 
compatible, apart from that of Kress (which is not to say that Kress quadrature is incompatible 
with fast summation; merely that a standard FMM will not work out of the box). 



Another important distinction is the one between (a) schemes in which the analytic split (1.2) 
into two smooth kernels must be explicitly known (i.e. the functions ip and ip are independently 
evaluated), and (b) schemes which merely need to access the overall kernel function k. The 
latter schemes are more flexible, since in applications the split is not always available (as in 



the axisymmetric example of section 7.4). However, as we will see, this flexibility comes with a 
penalty in terms of accuracy. 

The following schemes will be described: 

• Kapur— Rokhlin (section [3]) . This is the simplest scheme to implement, based upon an 
underlying periodic trapezoid rule. The weights, but not the node locations, are modified 
near the diagonal. No explicit split is needed. 

• Alpert (section |4]) . Also based upon the periodic trapezoid rule, and also not needing 
an explicit split, this scheme replaces the equi-spaced nodes near the diagonal with an 
optimal set of auxiliary nodes, at which new kernel evaluations are needed. 

• Modified Gaussian (section [5]) . Here the underlying quadrature is Gauss-Legendre 
panels, and new kernel evaluations are needed at each set of auxiliary nodes chosen for 
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each target node in the panel. These auxihary nodes are chosen using the algorithm of 
Kolm-Rokhlin [16]. No explicit split is needed. 

Kress (section [6]). This scheme uses an explicit split to create a spectrally-accurate 
product quadrature based upon the periodic trapezoid rule nodes. All of the matrix 



elements differ from the standard form (1.7), thus the scheme is not FMM-compatible 



We include it ctS Si benchmark where possible. 

Table [l] classifies these schemes (and a couple of others), according to whether they have under- 
lying global (periodic trapezoid rule) or panel-based quadrature, whether the split into the two 
smooth functions need be explicitly known or not, and whether they are FMM-compatible. 

In section [7] we present numerical tests comparing the accuracy of these quadratures in ID, 
2D, and 3D axisymmetric settings. We also demonstrate that some schemes have negative effects 
on the convergence rate in an iterative setting. We compare the advantages of the schemes and 
draw some conclusions in section [H 

1.3. Related work and schemes not compared. The methods described in this paper rely 
on earlier work |17[ [U [T^ [TC] describing high-order quadrature rules for integrands with weakly 
singular kernels. It appears likely that these rules were designed in part to facilitate Nystrom 
discretization of BIEs, but, with the exception of Kress PTT], the original papers leave most 
details out. (Kress describes the Nystrom implementation but does not motivate the quadrature 
formula; hence we derive this in section |6j) Some later papers reference the use (e.g. [HdH]) of 
high order quadratures but provide few details. In particular, there appears to have been no 
formal comparison between the accuracy of different approaches. 

There are several schemes that we do not have the space to include in our comparison. One of 
the most promising is the recent scheme of Helsing for Laplace [11] and Helmholtz [12j problems, 
which is panel-based but uses an explicit split in the style of Kress, and thus needs no extra 
kernel evaluations. We also note the recent QBX scheme [15] (quadrature by expansion) which 
makes use of off-curve evaluations and local expansions of the PDE. 

2. A BRIEF REVIEW OF LAGRANGE INTERPOLATION 

This section reviews some well-known (see, e.g., [2., Sec 3.1]) facts about polynomial interpo- 
lation that will be used repeatedly in the text. 

For a given set of distinct nodes {xj}jLi and function values {yj}jLi, the Lagrange interpo- 
lation polynomial L{x) is the unique polynomial of degree no greater than — 1 that passes 
through the A^ points {{xj, yj)}jLi. It is given by 

N 



where 



N 



(2.1) Hx) = l[ 



1=1 



Xi 



While polynomial interpolation can in general be highly unstable, it is perfectly safe and accurate 
as long as the interpolation nodes are chosen well. For instance, for the case where the nodes 
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{xj}j^i are the nodes associated with standard Gaussian quadrature on an interval / = [0,6], 
it is known [21 Thm. 3.2] that for any / G C^(/) 

N 



where 



c={ sup \f(^\s)\]/m 

\selo,b] J 



3. NySTROM discretization using the KAPUR-ROKHLIN QUADRATURE RULE 

3.1. The Kapur Rokhlin correction to the trapezoid rule. Recall that the standard 
A^+l-point trapezoid rule that approximates the integral of a function g G C°°[0, T] is h[g{0)/2 + 
g{h) + • • • + g{T — h) + g{T)/2], where the node spacing is 

h = ^, 

N ' 

and that it has only 2nd-order accuracy [18^ Thm. 12.1]. The idea of Kapur-Rokhlin [14J is to 
modify this rule to make it high-order accurate, by changing a small number of weights near the 
interval ends, and by adding some extra equi-spaced evaluation nodes Xj = hj for — /c < j < 
and N < j < N + k, for some small integer k > 0, i.e. nodes lying just beyond the ends. They 
achieve this goal for functions g G C°°[0,T], but also for the case where one or more endpoint 
behavior is singular, as in 



(3.1) 



g{x) = (p{x)s{x) + ipix) , 



where (p{x),ip{x) G C°°[0,T] and s(x) G C(0,r) is a known function with integrable singularity 
at zero, or at T, or at both zero and T. Accessing extra nodes outside of the interval suppresses 
the rapid growth with order of the weight magnitudes that plagued previous corrected trapezoid 
rules. However, it is then maybe unsurprising that their results need the additional assumption 
that ip,ip £ C^l-hk^T + hk]. 



Since we are interested in methods for kernels of the form (1.2), we specialize to a periodic 
integrand with the logarithmic singularity at x = (and therefore also at x = T, making both 
endpoints singular), 

XTT 



(3.2) 



(p{x) log 



sm ■ 



T 



+ 'il^{x) 



The mth-order Kapur-Rokhlin rule T, 
right endpoints is. 



m 



^ which corrects for a log singularity at both left and 



(3.3) 



'(9) = h 



?=-m 



lig{lh)+g{h) + g{2h) + ---+g{T 



h)+ j.eg{T + eh) 
e=-m 



Notive that the left endpoint correction weights {7^} 
the right endpoint. The convergence theorems from [1 



_^ are also used (in reverse order) at 
then imply that, for any fixed T-periodic 
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(3.4) 



g{x)dx -T^^^{g) 



0{h" 



as N 



oo 



The apparent drawback that function values are needed at nodes ih for — m < £ < —1, which 
he outside the integration interval, is not a problem in our case of periodic functions, since by 



periodicity these function values are known. This enables us to rewrite (3.3) as 



(3.5) T^^Hg) = h 



Y,{ll + l-t)9{^h)+g{h)+g{2h) + - ■ ■+g{T-h)+ ^ {^^ + ^^^)g{T + £h) 



which involves only the — 1 nodes interior to [0,T]. 

For orders m = 2,6, 10 the values of 7^ are given in the left-hand column of |14t Table 6]. 
In our periodic case, only the values 7^ + 7_f are needed; for convenience we give them in 
appendix [A] Notice that they are of size 2 for m = 2, of size around 20 for m = 6, and of size 
around 400 for m = 10, and alternate in sign in each case. This highlights the statement of 
Kapur-Rokhlin that they were only partially able to suppress the growth of weights in the case 
of singular endpoints |14j . 



3.2. A Nystrom scheme. We now can construct numbers ajj such that (1.4) holds. We 
start with the underlying periodic trapezoid rule, with equi-spaced nodes {xj}^;^ with Xj = hj, 
h = T/N. We introduce a discrete offset function between two nodes Xj and xj defined 

by 

Hhj) = 3-i (mod iV), -N/2 < < N/2 , 

and note that for each i G {1, . . . , A''}, and each \i\ < N/2, there is a unique j E {1, . . . , A^} such 
that = I. Two nodes Xi and Xj are said to be "close" if |^(i,j)| < m. Moreover, we call 

Xi and Xj "well-separated" if they are not close. 



We now apply (3.5) to the integral (1.4), which by periodicity of a and the kernel, we may 



rewrite with limits Xi to + T so that the log singularity appears at the endpoints. We also 
extend the definition of the nodes Xj = hj for all integers j, and get. 



Jo 



k{xi, x')a{x') dx' 



k{xi, x')a{x') dx' 



Xi 



i+N-l 

h k{xi,Xj)a{xj) 
j=i+i 



h X] (7^ + ^-.e)k{xi,Xi+e) a{xi+i) 



Wrapping indices back into {1, . . . , N}, the entries of the coefficient matrix A are seen to be, 
(3.6) 



if i = j, 

if Xi and xj are "well-separated" , 
if Xi and xj are "close", and i ^ j- 





+ 7e{i,j) + H^i, Xj) 

Notice that this is the elementwise product of a circulant matrix with the kernel matrix k{xi,Xj), 
and that diagonal values are ignored. Only 0{N) elements (those closest to the diagonal) differ 
from the standard Nystrom formula ( |1.7[ ). 
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Figure 2. Example of Alpert quadrature scheme of order / = 10 on the interval 
[0,1]. The original trapezoid rule had 20 points including both endpoints, i.e. 

= 19 and /i = 1/19. Correction nodes {Xp^l^i and {1 — Xpf^}^=i for m = 10 
and a = 6, are denoted by stars. 



4. Nystrom discretization using the Alpert quadrature rule 
4.1. The Alpert correction to the trapezoid rule. Alpert quadrature is another correction 



to the trapezoid rule that is high-order accurate for integrands of the form (3.1) on (0,r). The 
main difference with Kapur-Rokhlin is that Alpert quadrature uses node locations Ojff the equi- 
spaced grid xj = hj, but within the interval (0,T). Specifically, for the function g in (3.2) with 
log singularities at both ends, we denote by S^^^{g) the Zth-order Alpert quadrature rule based 
on an A+l-point trapezoid grid, defined by the formula 



N-a 



(4.1) 



(g) 



h^Wpgixph) + g{jh) + h^Wpg{T ■ 

p=i 



Xph). 



p=i 



j=a 



There are A^ — 2a + 1 internal equi-spaced nodes with spacing h = T/N and equal weights; 
these are common to the trapezoid rule. There are also m new "correction" nodes at each end 
which replace the a original nodes at each end in the trapezoid rule. The label "/th-order" is 
actually slightly too strong: the scheme is proven p3, Cor. 3.8] to have error convergence of order 
0{h''\ log /i|) as /i — 7- 0. The number m of new nodes needed per end is either l—l or For each 
order I, the integer a is chosen by experiments to be the smallest integer leading to positive 
correction nodes and weights. The following table shows the values of m and a for log-singular 
kernels for convergence orders / = 2, 6, 10, 16: 



Convergence order 


/i^ log/i 


h^\logh\ 


/iiO|log/i| 


h^^\logh\ 


Number of correction points m 


1 


5 


10 


15 


Width of correction window a 


1 


3 


6 


10 



The corresponding node locations xii • • • Xm and weights wi, ... Wm are listed in Appendix [Bj 
and illustrated for the case / = 10 in Figure [2] Details on how to construct these numbers by 
solving a nonlinear system can be found in [U Sec. 5], and the quadratures listed are adapted 
from |1, Table 8]. 



4.2. A Nystrom scheme. Recall that we wish to construct a matrix Uij that when applied to 
the vector of values {a{xj)}jLi approximates the action of the integral operator on the function 
o", evaluated at each target node Xj. To do this via the /th-order Alpert scheme with parameter 



a, we start by using periodicity to shift the domain of the integral in (1.4) to {xi,Xi + T), as 



in section (3.2). Since the 2m auxiliary Alpert nodes lie symmetrically about the singularity 

and 



location Xi, for simplicity we will treat them as a single set by defining Xp+ 



-Xp 
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Wp+m = Wp iov p = 1, . . . ,m. Then the rule (4.1 ) gives 
(4.2) 

N—a 2m 

k{xi,x')a{x') dx h"^ k{xi,Xi + ph) a{xi + ph) + ^ ^ Wpk{xi, Xi + Xph) cr{xi + Xph) ■ 

p=a p=l 

The values {Xp}p=i are not integers, so no auxiliary nodes coincide with any equispaced nodes 
{^^jljLi a-t which the vector of a values is given. Hence we must interpolate a to the auxiliary 
nodes {xi + Xp^}p=i- We do this using local Lagrange interpolation through M equispaced 
nodes surrounding the auxiliary source point xt + Xph- For M > / the interpolation error is 
higher order than that of the Alpert scheme; we actually use M = / + 3 since the error is then 
negligible. 

Remark 4.1. While high-order Lagrange interpolation through equi-spaced points is generally a 
bad idea due to the Runge phenomenon [23] , here we will be careful to ensure that the evaluation 
point always lies nearly at the center of the interpolation grid, and there is no stability problem. 

For each auxiliary node p, let the nth interpolation node index offset relative to i be 

o(f) := lxp-M/2\+n, 

and let the whole set be O^p) := {o^f^}^£i. For q e 0^p\ let the function np{q) := q- [xp-M/2\ 
return the node number of an index offset of q. Finally, let 

k=l \0n - J 

be the nth Lagrange basis polynomial for auxiliary node p. Applying this interpolation in a 
gives for the second term in ( |4.2[ ), 

2m M 

h ^Wpk{xi,Xi + Xph) ^L^n\xp)(^{Xi^^{p)) ■ 

p=l n=l 

Note that all node indices in the above will be cyclically folded back into the set {!,..., A^}. 
Recalling the notation i{i,j) from section 3.2 we now find the coefficient matrix A has entries 

(4.3) Qij = bij + dj, 

where the first term in (4.2) gives 

'0 if |£(i,j)| <«, 



^^■^^ " 1 hk{xi,Xj) if \i{i,j)\ > a, 

the standard Nystrom matrix (1.7) with a diagonal band set to zero, and the auxiliary nodes 
give 

2m 

(4.5) Cij = h ^ Wpk{xi,Xi + Xph)L''^^(^^^^jy~^{xp) . 

p=i 

Notice that the bandwidth of matrix does not exceed a + M/2, which is roughly /, and thus 
only 0{N) elements of Uij differ from those of the standard Nystrom matrix. 
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5. NySTROM DISCRETIZATION USING MODIFIED GAUSSIAN QUADRATURE 

We now turn to a scheme with panel-based underlying quadrature, namely the composite 
Gauss-Legendre rule. We recall that for any interval / = [0,6], the single-panel n-point Gauss- 
Legendre rule has nodes {xj}"^^ C I and weights {wj}^^^ C (0,oo) such that the identity 

(5.1) / fix)dx = y2wjf{xj) 

Jo ,=1 

holds for every polynomial / of degree at most 2n — 1. For analytic / the rule is exponentially 
convergent in n with rate given by the largest ellipse with foci and b in which / is analytic 
Ch. 19]. 



5.1. Modified Gaussian quadratures of Kolm Rokhlin. Suppose that given an interval 
[0, b] and a target point t £ [0, 6], we seek to approximate integrals of the form 

(5.2) [\ipis)S{t,s)+i^{s))ds, 







where ip and ■0 are smooth functions over [0, b] and S{t, s) has a known singularity as s — )• t. 
Since the integrand is non-smooth, standard Gauss-Legendre quadrature would be inaccurate if 



applied to evaluate (5.2). Instead, we seek an m-node "modified Gaussian" quadrature rule with 



weights {'Vfc}^i and nodes {yk}]^=i C [0,6] that evaluates the integral (5.2) to high accuracy. 
In particular, we use a quadrature formula of the form 

pb ™ 
(5.3) / {ip{s)S{t,s) + i'{s))ds^^Vk{(p{yk)S{t,yk) + 'ipiyk)) 

k=i 

which holds when ip and ■0 are polynomials of degree n. It is crucial to note that {'WfcjfcLi 
{yk}k=i depend on the target location t; unless t values are very close then new sets are needed 
for each different t. 



Next consider the problem of evaluating (5.2) in the case where the target point t is near the 



region of integration but not actually inside it, for instance t S [—6,0) U (6, 26]. In this case, 
the integrand is smooth but has a nearby singularity: this reduces the convergence rate and 
means that its practical accuracy would be low for the fixed n value we prefer. In this case, we 



can approximate the integral (5.2) with another set of modified Gaussian quadrature weights 



{vk}k=i aii'i nodes {yk}k=i [0;^] giving a quadrature formula analogous to (5.3). In fact, such 
weights and nodes can be found that hold to high accuracy for all targets t in intervals of the 
form [-W-P+\-W-P]. 

Techniques for constructing such modified quadratures based upon nonlinear optimization are 
presented by Kolm-Rokhlin [16] . In the numerical experiments in section [7| we use a rule with 
n = 10, m = 20, and m' = 24. This rule leads to fairly high accuracy, but is not of so high order 
that clustering of the end points becomes an issue in double precision arithmetic. Since the full 
tables of quadrature weights get long in this case, we provide them as text files at jlOj. 

5.2. A Nystrom scheme. We partition the domain of integration as 

Np 

[o,T] = \Jnp, 
p=i 
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where the Op's are non-overlapping subintervals cahed panels. For simplicity, we for now assume 



that the panels are equi-sized so that ilp 



Note that an adaptive version would 



T(p-l) Tp 
Np ' Np 

have variable-sized panels. On each panel, we place the nodes of an n-point Gaussian quadrature 
to obtain a total of = nNp nodes. Let {xi}^^ and {'Wi}^^ denote the nodes and weights of 
the resulting composite Gaussian rule. 

Now consider the task of approximating the integral (1.4) at a target point Xj. We decompose 
the integral as 



(5.4) 



Np 

') cj(x') fix' = / k{xi,x') a{x') dx . 

q=l 



We will construct an approximation for each panel-wise integral 



(5.5) 



/c(xj, x)a{x) dx' , 



expressed in terms of the values of a at the Gaussian nodes in the source panel Qq. There are 
three distinct cases: 

Case 1: Xi belongs to the source panel Vtq: The integrand is now singular in the domain of 
integration, but we can exploit that a is still smooth and can be approximated via polynomial 
interpolation on this single panel. Using the set of n interpolation nodes {x^ ■ x^ G ^q}, let Lj 
be the Lagrange basis function corresponding to node j. Then, 



(5.6) 



aix') 



Yl Lj{x')a{xj) . 



Inserting (5.6) into the integral in (5.5) we find 



k{xi, x')(j{x') dx' 




k{xi,x')Lj{x') dx' \ a{xj). 



Let {vi^k}^=i and {yi,k}^=i be the modified Gaussian weights and nodes in the rule (5.3) on the 
interval Vtq associated with the target t = Xi. Using this rule, 

(5.7) / k{x^,x')a{x')dx' ^ ^ 'yZvi,kk{xi,yi^k) Lj{yi^k) cr{xj). 

■^^1 j : Xj^q \k=l / 

Note that the expression in brackets gives the matrix element a^j-. Here the auxiliary nodes yj^fc 
play a similar role to the auxiliary Alpert nodes Xi + Xph from section |4j new kernel evaluations 
are needed at each of these nodes. 

Case 2: Xi belongs to a panel fip adjacent to source panel 0,q: In this case, the kernel k is smooth, 
but has a singularity closer to Qq than the size of one panel, so standard Gaussian quadrature 
would still be inaccurate. We therefore proceed as in Case 1: we replace a by its polynomial 
interpolant and then integrate using the modified quadratures described in Section 5.1 , The end 
result is a formula similar to (5.7) but with the sum including m' rather than m terms, and with 
Vi^k and iji k replacing Vi^k and yi^k, respectively. 
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Case 3: Xi is well-separated from the source panel 0,q: By "well-separated" we mean that Xi and 
0,q are at least one panel size apart in the parameter x. (Note that if the curve geometry involves 
close-to-touching parts, then this might not be a sufficient criterion for being well-separated in 
M?; in practice this would best be handled by adaptivity.) In this case, both the kernel k and 
the potential a are smooth, so the original Gaussian rule will be accurate. 




Combining (5.8) and (5.7), we find that the Nystrom matrix elements a^.j are given by 

XlfcLi Vi^k k{xi,yi^k) Lj{yi^k), if Xi and xj are in the same panel, 
YlT=i ^i,k K^i,yi,k) Ljivi^k), if Xi and xj are in adjacent panels, 
k{xi, Xj) Wj, if Xi and Xj are in well-separated panels. 



6. Nystrom discretization using the Kress quadrature rule 
The final scheme that we present returns to an underlying periodic trapezoid rule, but de- 



mands separate knowledge of the smooth kernel functions ip and appearing in (1.2). We first 
review spectrally-accurate product quadratures, which is an old idea, but which we do not find 
well explained in the standard literature. 

6.1. Product quadratures. For simplicity, and to match the notation of \I7j, we fix the period 
T = 2tt and take to be even. The nodes are thus Xi = 2TTi/N , i = 1, . . . , N . 

A product quadrature approximates the integral of the product of a general smooth 27r-periodic 
real function / with a fixed known (and possibly singular) 27r-periodic real function g, by a 
periodic trapezoid rule with modified weights Wj, 

(6.1) / fis)gis)ds ^ ^w,f{x,) . 







Using the Fourier series /(s) = ^^g^ fn^^^^, and similar for g, we recognize the integral as an 
inner product and use Parseval, 

f-2n 

(6.2) / fis)gis)ds = 271^2 fng^ 

Since / is smooth, decays to zero with a high order as \n\ — oo. Thus we can make two 
approximations. Firstly, we truncate the infinite sum to Yl'\n\<N/2^ where the prime indicates 
that the extreme terms n = ±iV/2 are given a factor of 1/2. Secondly, we use the periodic 
trapezoid rule to evaluate the Fourier coefficents of /, i.e. 

r2n 1 N 



(6.3) = ^1 '^^"^"V(^)ds « ^Ee"^'^'V(x,) . 

Although the latter introduces aliasing (one may check that the latter sum is exactly fn + fn+ N + 
fn-N + fn+2N + fn-2N + ' ' ' )) ^hc decay of fn means that errors decay to high order with N. 
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Substituting (6.3) into the truncated version of (6.2) gives 
(6.4) / f{s)g{s)ds ^2^ 5^]^ E ^"'"'^ /(^^•) ■ 

\n\<N/2 j=l 



.7 = 1 ^ 



\n\<N/2 



The bracketed expression gives the weights in (6.1). Since g is real (hence g-n = g-n) 
(6.5) 

N/2-1 



W 



2tt sr-^ 



H<N/2 



27r 
iV 



<7o+ E 2Re(ff„e--0 + Re(<7^/2e^^-^/2) 



n=l 



j = l,...,N. 



6.2. The Kress quadrature. To derive the scheme of Kress (originaUy due to Martensen- 
Kussmaul; see references in [17j) we note the Fourier series (proved in |181 Thm. 8.21]), 



(6.6) 



g{s) = log (4 sin 



0, n = 0, 

-l/\n\, 



Translating g by a displacement t S M corresponds to multiplication of gn by e * . Substituting 
this displaced series into (6.5) and simplifying gives 

r.27r 

(6.7) 



p2TT / + \ ^ 



'3 1 ' 



where the weights, which depend on the target location t, are 



(6.8) Rf'^it) 



'iV 



.N/2-l 



^ 1 , ^ I N , 

2^ - cos n{xj -t) + cos —{xj - t) 



n=l 



j = l,...,N. 



This matches pTl (3.1)] (note the convention on the number of nodes differs by a factor 2). 
When a smooth function is also present, we use the periodic trapezoid rule for it, to get 



(6.9) / log Usin^^ ^{s)+^l.{s)ds ^ ^ i?f/2)(t) ^(x,) + ^ ^ ^(x,) . 

Assuming the separation into <p and tp is known, this gives a high-order accurate quadrature; in 
fact for ip and ip analytic, it is exponentially convergent |17| . 

6.3. A Nystrom scheme. We use the above Kress quadrature to approximate the integral 
(1.4) where the kernel has the form (1.2) with T = 2tt, and the functions (p{x,x') and 'iIj{x,x') 
are separately known. Applying (6.9), with h = 2tt/N, gives 

r2n N N 

(6.10) / k{xi,x')a{x') dx' w ^.^j {xi)^{xiiXj)a{xj) + hy ^'ij)[xi,Xj)a{xj) . 

Using the symbol Rj^^"^^ := R^'^^'^\o), and noticing that Rj^^'^\xi) depends only on \i — j\, we 
find that the entries of the coefficient matrix A are. 



N 



(6.11) 



i,j = R^^fJ^^^^ix^Xj) + h'tp{xi,Xj) 
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Note that R 



(N/2) 



is a dense circulant matrix, and all A'^^ elements differ from the standard 



Nystrom matrix (1.7). Since and ^p do not usually have fast potential-theory based algorithms 



to apply them, the Kress scheme is not FMM-compatible. 



7. Numerical experiments 

We now describe numerical experiments in ID, 2D and 3D applications that illustrate the 
performance of the quadratures described in sections 3][6 The experiments were carried out on 
a Macbook Pro with 2.4GHz Intel Core 2 Duo and 4GB of RAM, and executed in a MATLAB 



environment. Once the Nystrom matrix is filled, the linear system (1.6) is solved via MATLAB's 



backslash (mldivide) command. In all examples below, the errors reported are relative errors 
measured in the L°° -norm, \\u^ — ^||oo/||^||oo; where u is the reference solution and is the nu- 
merical solution. For each experiment we compare the performance of the different quadratures. 
Specifically, we compare the rules of Kapur-Rokhlin of orders 2, 6, and 10; Alpert of orders 2, 
6, 10, and 16; modified Gaussian with n = 10 points per panel; and (where convenient) Kress. 
Our implementation of the modified Gaussian rule uses m = 20 auxiliary nodes for source and 
target on the same panel, and m' = 24 when on adjacent panels. 

The quadrature nodes and weights used are provided in appendices and at the website pCU] . 



10" 



10"' 



10 



10 



10 



10" 



10" 



10" 




■ mod. gauss. 



2"" K-R 



10* K-R 
- 2""^ Alpert 
6* Alpert 
10* Alpert 
16* Alpert 



Kress 



10^ 

N (number of nodes) 



10" 



Figure 3. Error results for solving the integral equation (7.1) in Section 7.1 



7.1. A ID integral equation example. We solve the one-dimensional integral equation 
(7.1) u{x)+ I k{x,x)u{x')dx = f{x), x £ [0, 2tt] 
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associated with a simple kernel function having a periodic log singularity at the diagonal, 

X — a 



(7.2) k{x,x') = ^log 



sm ■ 



-logUsm^— j - -log2 



thus the smooth functions (/7(x, x') = 1/4 andijj{x,x') = —(1/2) log2 are constant. This kernel is 



similar to that arising from the Laplace single- layer operator on the unit circle. Using (6.6) one 
may check that the above integral operator has exact eigenvalues — 7rlog2 (simple) and — 7r/(2n), 
n = 1,2, .. . (each doubly-degenerate). Thus the exact condition number of the problem ( |7.1[ ) 
is ((7rlog2) — 1)/(1 — ~ 5.5. We choose the real-analytic periodic right-hand side f{x) = 
sin(3a;) e'^°''(^^). The solution u has ||it||oo ~ 6.1. We estimate errors by comparing to the Kress 



solution at = 2560. (In passing we note that the exact solution to (7.1) could be written 
analytically as a Fourier series since the Fourier series of / is known in terms of modified Bessel 
functions.) When a solution is computed on panel-based nodes, we use interpolation back to 
the uniform trapezoid grid by evaluating the Lagrange basis on the n = 10 nodes on each panel. 

In Figure [7} the errors in the L°°-norm divided by Halloo are presented for N = 20, 40, 80, . . . , 1280. 
We see that the rules of order 2, 6, and 10 have the expected convergence rates, but that Alpert 
has prefactors a factor 10^ to 10^ smaller the Kapur-Rokhlin. We also see that Kress is the most 
efficient at any desired accuracy, followed by the three highest-order Alpert schemes. These four 
schemes flatten out at 13 digits, but errors start to grow again for larger A^, believed due to the 
larger linear system. Note that modified Gaussian performs roughly as well as 6th-order Alpert 
with twice the number of points, and that it flattens out at around 11 digits. 

7.2. Combined field discretization of the Helmholtz equation in M?. In this section, we 
solve the Dirichlet problem for the Helmholtz equation exterior to a smooth domain Q C 
with boundary F, 

(7.3) -Ait-w\ = 0, in E = n'', 

(7.4) u = f, on F, 

where w > is the wavenumber. u satisfies the Sommerfeld radiation condition 

(7.5) lim r^/^ -iuju] =0, 



where r = and the limit holds uniformly in all directions. A common approach [7", Ch. 3] is 
to represent the solution to (|7.3|) via both the single and double layer acoustic potentials. 



u{x) = j k{x,x') a{x') dl{x') 



(7.6) = [^^^^^ - H^, ^')) ^(^0 dl{x'), XGE, 

where (j){x,x') = ^Hq'^\uj\x — x'\) and Hq^^ is the Hankel function of the first kind of order 
zero; n is the normal vector pointing outward to F, and dl the arclength measure on F. The 



motivation for the combined representation (7.6) is to obtain the unique solvability to problem 
( |7.3 7.4) for all oj > 0. The corresponding boundary integral equation we need to solve is 



(7.7) ]-a[x)+ f k{x,x')a{x')dl{x') = f{x), a; e F, 

2 Jy 
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(a) 0.5 A diameter 
(w = 2.8) 



(b) 5 A diameter 
(w = 28) 



(c) 50 A diameter 
(w = 280) 



- mod. gauss. 



- 'i^ - 10 K-R 
— ' — 2"" Alpert 
- 6* Alpert 



10* Alpert 
16* Alpert 
Kress 




600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 
N (number of nodes) 



Figure 4. Error results for the exterior planar Helmholtz problem (7.3) in Sec- 
tion 7.2 solved on the starfish domain of Figure [T} 
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where k{x, x') = ^^^j^ - iuj (pix, x') 



To convert (7.7) into an integral equation on the real line, we need to parametrize F by a 
vector-valued smooth function r : [0,T] — t- M^. By changing variable, (7.7) becomes 

(7.8) ^aiT{t)) + £ k{T{t),Tis))aiTis))\dT/ds\ds = fiT{t)), te[0,T]. 

To keep our formula uncluttered, we rewrite the kernel as 

(7.9) m{t,s) = k{T{t),T{s)) \dT/ds\, 
as well as the functions 

a{t) = a{r{t)) and /(t) = /(T(t)). 
Thus we may write the integral equation in standard form, 

(7.10) a{t) + 2 [ m{t,s)a{s)ds = 2f{s), t £ [0,T], 

Jo 

and apply the techniques of this paper to it. 

Remark 7. 1 . All the quadrature schemes apart from that of Kress are now easy to implement by 
evaluation of m(t,s). However, to implement the Kress scheme, separation of the parametrized 
single- and double-layer Helmholtz kernels into analytic functions if and ip is necessary, and not 
trivial. We refer the reader to [TTl Sec. 2] or [7, Ch. 3]. 

We assess the accuracy of each quadrature rule for the smooth domain shown in Figure , 
varying N and the wavenumber uj. Specifically, we varied wavenumbers such that there are 0.5, 
5 and 50 wavelengths across the domain's diameter. The right-hand side is generated by a sum 
of five point sources inside il., with random strengths; thus the exact exterior solution is known. 
Errors are taken to be the maximum relative error at the set of measurement points shown in 
Figure [l][|c). Notice that sources and measurement points are both far from F, thus no challenges 
due to close evaluation arise here. 

The results are shown in Figure |4] At the lowest frequency, results are similar to Figure [7j 
except that errors bottom out at slightly higher accuracies, probably because of the smoothing 
effect of evaluation of u at distant measurement points. In terms of the error level each scheme 
saturates at, Kress has a 2-3 digit advantage over the others at all frequencies. At the highest 
frequency, Figure Qc), there are about 165 wavelengths around the perimeter F, hence the 
N = 1000 at which Kress is fully converged corresponds to around 6 points per wavelength. 
At 10 points per wavelength (apparently a standard choice in the engineering community, and 
roughly halfway along the N axis in our plot), 16th-order Alpert has converged at 10 digits, 
while modified Gaussian and lOth-order Alpert are similar at 8 digits. The other schemes are 
not competitive. 

7.3. Effect of quadrature scheme on iterative solution efficiency. When the number 
of unknowns N becomes large, iterative solution becomes an important tool. One standard 
iterative method for nonsymmetric systems is GMRES [22j. In Table [2] we compare the numbers 



of GMRES iterations needed to solve the linear system (1.4) arising from the low- frequency 



Helmholtz EVP in section 7.2 (0.5 wavelengths across), when the matrix was constructed via 
the various quadrature schemes. We also give the condition number of the system. We see that 
most schemes result in 14 iterations and a condition number of around 3.5; this reflects the 
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scheme 


mod. Gauss 


2nd K-R 


6th K-R 


10th K-R 


2nd Alpert 


6th Alpert 


10th Alpert 


16th Alpert 


Kress 


cond # 


3.95 


3.52 


3.68 


169 


3.52 


3.52 


3.52 


3.52 


3.52 


# iters 


14 


14 


22 


206 


14 


14 


14 


14 


14 



Table 2 . Condition numbers of the Nystrom system matrix ( ^ I + A) , and num- 
bers of GMRES iterations to reach residual error 10~^^, for ah the quadrature 
schemes. 



(^) Imodified bauss 
2"^ K-R 
6* K-R 
10'" K-R 
2™ Alpert 
6* Alpert 
10 ^ Alpert 
16'" Alpert 
Kress 
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-0.6 
-0.8 



0.2 0.4 0.6 0.8 1 1.2 



Figure 5. (a) Magnitude of eigenvalues of the matrix (^l + A) associated 
with the Nystrom discretization of the Helmholtz BVP (7.7). The system size 
is = 640 and the wave- number cj corresponds to a contour of size 0.5 wave- 
lengths, (b) Eigenvalues in the complex plane associated with lOth-order Kapur- 
Rokhlin (dots) and Kress (crosses) quadratures, (c) Same plot as (b), but zoomed 
in to the origin. 



fact that the underlying integral equation is Fredholm 2nd-kind and well-conditioned. However, 
6th-order and particularly lOth-order Kapur-Rokhlin require many more iterations (a factor 
15 more in the latter case), and have correspondingly high condition numbers. In a practical 
setting this would mean a much longer solution time. 

In order to understand why, in Figure [5] we study the spectra of the system matrix; since the 
operator is of the form Id/2 -|- compact, eigenvalues should cluster near 1/2. Indeed this is the 
case for all the schemes. However, 6th- and lOth-order Kapur-Rokhlin also contain additional 
eigenvalues that do not appear to relate to those of the operator. In the lOth-order case, these 
create a wide "spray" of many eigenvalues, which we also plot in the complex plane in (b) and 
(c). A large number (around 200) of eigenvalues fall at much larger distances than the true 
spectrum, and on both sides of the origin; we believe they cause the slow GMRES convergence. 
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The corresponding eigenvectors are oscillatory, typically alternating in sign from node to node. 
We believe this pollution of the spectrum arises from the large alternating weights 7; in these 
schemes. Note that one spurious eigenvalue falls near the origin; this mechanism could induce an 
arbitrarily large condition number even though the integral equation condition number is small. 
Although we do not test the very large values where iterative methods become essential, we 
expect that our conclusions apply also to larger N. 





(a) 





(b) 



Figure 6. Domains used in numerical examples in Section 7.4 All items are 
rotated about the vertical axis, (a) A sphere, (b) A starfish torus. 



7.4. The Laplace BVP on axisymmetric surfaces in R^. In this section, we compare 
quadratures rules applied on kernels associated with BIEs on rotationally symmetric surfaces in 
M^. Specifically, we considered integral equations of the form 

(7.11) a{x) + j k{x,x')a{x')dA{x') = f{x), xeT, 

under the assumptions that F is a surface in obtained by rotating a curve 7 about an axis 
and the kernel function k is invariant under rotation about the symmetry axis. Figure [6] depicts 
domains used in numerical examples: the generating curves 7 are shown in the left figures 



and the axisymmetric surfaces F are shown in the right ones. The BIE (7.11) on rotationally 
symmetric surfaces can via a Fourier transform be recast as a sequence of equations defined on 
the generating curve in cylindrical coordinates, i.e. 



(7.12) cr„(r, z) + V27r / kn{r, z,r' , z') an{r' , z') r' dl{r' , z') = fn{r, z), {r,z)ej, n E Z, 



r 
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(a) 10° 




N (number of nodes) 



(b) 10° 




N (number of nodes) 

Figure 7. Error results for the 3D interior Dirichlet Laplace problem from see 
tion 



7.4 solved on the axisymmetric domains (a) and (b) respectively shown in 
Figure [6j 
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where an, fn, and kn denote the Fourier coefficients of a, f, and k, respectively. Details on 
how to truncate the Fourier series and construct the coefficient matrices for Laplace problem 
and Helmholtz problem can be found in [2B]. In the following experiments, we consider the BIE 
( |7.11 ) which arises from the interior Dirichlet Laplace problem, in which case 

n{x') ■ {x — x') 



(7.13) k{x,x') 



AttIx — x'P 



As we recast the BIE defined on P to a sequence of equations defined on the generating curve 
7, it is easy to see that the kernel function kn has a logarithmic singularity as {r',z') — )• (r, z). 
In this experiment, 101 Fourier modes were used. We tested all the quadrature schemes apart 
from that of Kress, since we do not know of an analytic split of the axisymmetric kernel kn into 
smooth parts ip and ip. 



Equation (7.11 ) was solved for Dirichlet data / corresponding to an exact solution u generated 
by point charges placed outside the domain. The errors reported reflect the maximum of the 
point-wise errors (compared to the known analytic solution) sampled at a set of target points 
inside the domain. 

The results are presented in Figure [Tj The most apparent feature in (a) is that, because 
the curve 7 is open, the schemes based on the periodic trapezoid rule fail to give high-order 
convergence; rather, it appears to be approximately 3rd-order. Panel-based schemes are able 
to handle open intervals as easily as periodic ones, thus modified Gaussian performs well: it 
reaches 12-digit accuracy with only around = 100 points. (We remark that the problems 
associated with an open curve are in this case artificial and can be overcome by a better problem 
formulation. We deliberately chose to use a simplistic formulation to simulate an "open curve" 
problem.) In (b), all functions are again periodic since 7 is closed; modified Gaussian performs 
similarly to the three highest-order Alpert schemes with around 1.5 to 2 times the number of 
points. 



8. Concluding remarks 

To conclude, we make some informal remarks on the relative advantages and disadvantages 
of the different quadrature rules that we have discussed. Our remarks are informed primarily 
by the numerical experiments in section [Tj 

Comparing the three schemes based upon nodes equi-spaced in parameter (Kapur-Rokhlin, 
Alpert, and Kress), we see that Kress always excels due to its superalgebraic convergence, 
converging fully at around 6 points per wavelength at high frequency, and its small saturation 
error of 10~^^ to 10~^^. However, the analytic split required for Kress is not always available 
or convenient, and Kress is not amendable to standard FMM-style fast matrix algebra. Kapur- 
Rokhlin and Alpert show their expected algebraic convergence rates at orders 2, 6, and 10, and 
both seem to saturate at around 10^^^. However, Alpert outperforms Kapur-Rokhlin since its 
prefactor is much lower, resulting in 2-8 extra digits of accuracy at the same A^. Another way to 
compare these two is that Kapur-Rokhlin requires around 6 to 10 times the number of unknowns 
as Alpert to reach comparable accuracy. The difference in a high-frequency problem is striking, 
as in Figure [4| The performance of 16th-order Alpert is usually less than 1 digit better than 
lOth-order Alpert, apart from at high frequency when it can be up to 3 digits better. 

Turning to the panel-based modified Gaussian scheme, we see that in the low-frequency set- 
tings it behaves like lOth-order Alpert but requires around 1.5 to 2 times the to reach similar 
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accuracy. This may be related to the fact that Gauss-Legendre panels would need a factor 7r/2 
higher than the trapezoid rule to achieve the same largest spacing between nodes; this is 
the price to pay for a panel-based scheme. However, for medium and high frequencies, 10th- 
order Alpert has little distinguishable advantage over modified Gaussian. Both are able to reach 
around 10 digit accuracy at 15 points per wavelength. Modified Gaussian errors seem to saturate 
at around 10~^^ to 10^^^. It therefore provides a good all-round choice, expecially when adap- 
tivity is anticipated, or global parametrizations are not readily constructed. One disadvantage 
relative to the other schemes is that the auxiliary nodes require kernel evaluations that are very 
close to the singularity (10^^ or less; for Alpert the minimum is only around 10^^). 

We have not tested situations in which adaptive quadrature becomes essential; in such cases 
modified Gaussian would excel. However, a hint of the convenience of modified Gaussian is given 
by its effortless handling of an open curve in Figure [7]^a) where the other (periodic) schemes 
become low-order (corner-style reparametrizations would be needed to fix this [7, Sec. 3.5]). 

In addition, we have showed that, in an iterative solution setting, higher-order Kapur-Rokhlin 
can lead to much slower GMRES convergence than any of the other schemes. We believe this 
is because it introduces many large eigenvalues into the spectrum, unrelated to those of the 
underlying operator. Thus lOth-order Kapur-Rokhlin should be used with caution. With that 
said, Kapur-Rokhlin is without doubt the simplest to implement of the four schemes, since no 
interpolation or new kernel evaluations are needed. 

We have chosen to not report computational times in this note since our MATLAB imple- 
mentations are far from optimized for speed. However, it should be mentioned that both the 
Alpert method and the method based on modified Gaussian quadrature require a substantial 
number (between 20A^ and SON in the examples reported) of additional kernel evaluations. 

For simplicity, in this note we limited our attention to the case of smooth contours, but both 
the Alpert and the modified Gaussian rule can with certain modifications be applied to contours 
with corners, see, e.g., |2T1 [5l El SJ [l3l [l2]. We plan to include the other recent schemes shown 
in Table [T| and curves with corners, in future comparisons. 
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Appendix A. Tables of Kapur-Rokhlin Quadrature weights 



2"''-order Kapur-Rokhlin correction weights 
integrals of the form f f{x) + g{x) log(a;) dx 


INDEX I 


WEIGHTS 7i + 7_i 


1 
2 


1.825748064736159e-F00 
-1.325748064736159e-F00 
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6 -order Kapur-Rokhlin correction weights 
integrals of the form f{x) + g{x) log(a;) dx 



INDEX I 



WEIGHTS 7( + 7- 



4.967362978287758e+00 

-1.620501504859126e+01 
2.585153761832639O+01 
-2.222599466791883C+01 
9.930104998037539e+00 
-1.817995878141594e+00 



10*'''-or(lor 


Kapiir Rokhhii comx:tioii weights 


iiit(;;j,ials ( 


f tlic form /jj^ /(•'■) — !■](■'') l^'Si-'') '^-'i' 


INDEX / 


WEIGHTS 7; + j-i 


1 


7.832432020568779e+00 


2 


-4.565161670374749e+01 


3 


1.452168846354677e+02 


4 


-2.901348302886379e+02 


5 


3.870862162579900e+02 


6 


-3.523821383570681e+02 


7 


2.172421547519342e+02 


8 


-8.707796087382991e+01 


9 


2.0535842660726356+01 


10 


-2. 166984103403823e+00 



Appendix B. Tables of Alpert Quadrature rules 



2"'°'-order Alpert Quadrature Rule for 
integrals of the form f{x) + g{x) log(x) dz, 
with a = 1 


NODES 


WEIGHTS 


1.591549430918953e-01 


5.000000000000000e-01 



6*'*-order Alpert Quadrature Rule for 
integrals of the form f{x) + g{x) log(2;) dx, 
with a = 3 


NODES 


WEIGHTS 


4.004884194926570e-03 

7.745655373336686C-02 
3.972849993523248e-01 
1.075673352915104e+00 
2.0037969271 11872e+00 


1.671879691147102e-02 

1.6369583714473606-01 
4.981856569770637C-01 
8.372266245578912e+00 
9.841730844088381e+00 
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10*'*-order Alpert Quadrature Rule for 
integrals of the form f{x) + g{x) log(x) dx, 
with a = 6 


NODES 


WEIGHTS 


1.175089381227308e-03 

1.877034129831289e-02 

9.686468391426860e-02 

3.004818668002884e-01 

6.901331557173356e-01 

1.293695738083659e+00 

2.090187729798780C+00 

3.016719313149212e+00 

4.001369747872486e+00 

5.000025661793423e+00 


4.560746882084207e-03 
3.810606322384757e-02 
1.293864997289512e-01 
2.884360381408835e-01 
4.958111914344961e-01 
7.077154600594529e-01 
8.741924365285083e-01 
9.661361986515218e-01 
9.957887866078700e-01 
9.998665787423845e-01 
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